Geometry and mesh for motor

Decide whether to draw the whole motor or just an 8th slice

pizza_slice = True

Define radii and points that help define the magnets, coils, and airgaps.

origin = (0,0)
# inner radius rotor
r1 = 26.5*10**(-3)
# outer radius rotor
r2 = 78.63225*10**(-3)
# sliding mesh rotor
r4 = 78.8354999*10**(-3)
# sliding mesh stator
r6 = 79.03874999*10**(-3)
# inner radius stator
r7 = 79.242*10**(-3)
# outer radius stator
r8 = 116*10**(-3)

# Points for magnet1 and air around magnet1
m1 = (69.23112999*10**(-3),7.535512*10**(-3))
m2 = (74.828958945*10**(-3),10.830092744*10**(-3))
m3 = (66.13621099700001*10**(-3),25.599935335*10**(-3))
m4 = (60.53713*10**(-3),22.30748*10**(-3))
a5 = (69.75636*10**(-3),5.749913*10**(-3))
a6 = (75.06735*10**(-3),3.810523*10**(-3))
a7 = (65.6868747*10**(-3),26.3184618*10**(-3))
# a7 = (65.3506200*10**(-3),26.51379*10**(-3))
a8 = (59.942145092*10**(-3),24.083661604*10**(-3))

# Points for magnet2 and air around magnet2
m5 = (58.579985516*10**(-3), 27.032444757*10**(-3))
m6 = (64.867251151*10**(-3),28.663475405*10**(-3))
m7 = (60.570096319*10**(-3),45.254032279*10**(-3))
m8 = (54.282213127*10**(-3),43.625389857*10**(-3))
a1 = (53.39099766*10**(-3),45.259392713*10**(-3))
a2 = (55.775078884*10**(-3),50.386185578*10**(-3))
a3 = (59.41521771*10**(-3),25.355776837*10**(-3))
a4 = (65.12210917100001*10**(-3),27.707477175*10**(-3))

# Points for Stator Nut and air in the stator
s1 = (79.04329892000*10**(-3),3.9538335974*10**(-3))
s2 = (80.143057128*10**(-3),4.0037794254*10**(-3))
s3 = (80.387321219*10**(-3),2.965459706*10**(-3))
s4 = (98.78501315600001*10**(-3),3.9007973292*10**(-3))
s5 = (98.44904989600001*10**(-3),9.026606148400001*10**(-3))
s6 = (80.086666706*10**(-3),7.5525611543*10**(-3))
s7 = (79.980020247*10**(-3),6.4912415424*10**(-3))
s8 = (78.88229587*10**(-3),6.4102654448*10**(-3))

Next, use the points to define drawing functions for the magnets and coils

import numpy as np
import netgen.occ as occ
from numpy import sin, cos, pi

def rotate(m,k,p):
    a = np.pi/p
    mx = m[0]*np.cos(k*a) -m[1]*np.sin(k*a)
    my = m[0]*np.sin(k*a) +m[1]*np.cos(k*a)
    return (mx, my, 0)

def drawMagnet1(k):
    m1new = rotate(m1,k,4); m2new = rotate(m2,k,4)
    m3new = rotate(m3,k,4); m4new = rotate(m4,k,4)

    a5new = rotate(a5,k,4); a6new = rotate(a6,k,4)
    a7new = rotate(a7,k,4); a8new = rotate(a8,k,4)

    #Draw magnet
    seg1 = occ.Segment(m1new,m2new); seg2 = occ.Segment(m2new,m3new)
    seg3 = occ.Segment(m3new,m4new); seg4 = occ.Segment(m4new,m1new)
    magnet1 = occ.Face(occ.Wire([seg1,seg2,seg3,seg4]))

    #Draw air around magnet
    air_seg1 = occ.Segment(m1new,a5new); air_seg2 = occ.Segment(a5new,a6new)
    air_seg3 = occ.Segment(a6new,m2new); air_seg4 = occ.Segment(m2new,m1new)
    air_magnet1_1 = occ.Face(occ.Wire([air_seg1,air_seg2,air_seg3,air_seg4]))

    air_seg5 = occ.Segment(m4new,m3new); air_seg6 = occ.Segment(m3new,a7new)
    air_seg7 = occ.Segment(a7new,a8new); air_seg8 = occ.Segment(a8new,m4new)
    air_magnet1_2 = occ.Face(occ.Wire([air_seg5,air_seg6,air_seg7,air_seg8]))

    return (magnet1,air_magnet1_1,air_magnet1_2)

def drawMagnet2(k):
    m5new = rotate(m5,k,4); m6new = rotate(m6,k,4)
    m7new = rotate(m7,k,4); m8new = rotate(m8,k,4)

    a1new = rotate(a1,k,4); a2new = rotate(a2,k,4)
    a3new = rotate(a3,k,4); a4new = rotate(a4,k,4)

    #Draw magnet
    seg1 = occ.Segment(m5new,m6new); seg2 = occ.Segment(m6new,m7new)
    seg3 = occ.Segment(m7new,m8new); seg4 = occ.Segment(m8new,m5new)
    magnet2 = occ.Face(occ.Wire([seg1,seg2,seg3,seg4]))

    #Draw air around magnet
    air_seg1 = occ.Segment(m5new,a3new); air_seg2 = occ.Segment(a3new,a4new)
    air_seg3 = occ.Segment(a4new,m6new); air_seg4 = occ.Segment(m6new,m5new)
    air_magnet2_1 = occ.Face(occ.Wire([air_seg1,air_seg2,air_seg3,air_seg4]))

    air_seg5 = occ.Segment(m8new,m7new); air_seg6 = occ.Segment(m7new,a2new)
    air_seg7 = occ.Segment(a2new,a1new); air_seg8 = occ.Segment(a1new,m8new)
    air_magnet2_2 = occ.Face(occ.Wire([air_seg5,air_seg6,air_seg7,air_seg8]))

    return (magnet2,air_magnet2_1,air_magnet2_2)

def drawStatorNut(k):

    s1new = rotate(s1,k,24); s2new = rotate(s2,k,24)
    s3new = rotate(s3,k,24); s4new = rotate(s4,k,24)
    s5new = rotate(s5,k,24); s6new = rotate(s6,k,24)
    s7new = rotate(s7,k,24); s8new = rotate(s8,k,24)

    #Draw stator coil
    seg1 = occ.Segment(s2new,s3new); seg2 = occ.Segment(s3new,s4new)
    seg3 = occ.Segment(s4new,s5new); seg4 = occ.Segment(s5new,s6new)
    seg5 = occ.Segment(s6new,s7new); seg6 = occ.Segment(s7new,s2new)
    stator_coil = occ.Face(occ.Wire([seg1,seg2,seg3,seg4,seg5,seg6]))

    #Draw air nut in the stator
    air_seg1 = occ.Segment(s1new,s2new); air_seg2 = occ.Segment(s2new,s7new)
    air_seg3 = occ.Segment(s7new,s8new); air_seg4 = occ.Segment(s8new,s1new)
    stator_air = occ.Face(occ.Wire([air_seg1,air_seg2,air_seg3,air_seg4]))

    stator_air = stator_air-(stator_air*air_gap_stator)
    return (stator_coil,stator_air)

Use the functions to define the shapes using the OCC technology.

domains = []

h_max = 1

h_air_gap = r6-r4 #0.05*h_max
h_air_magnets = h_max
h_coils = h_max
h_stator_air = h_max
h_magnets = h_max
h_stator_iron = h_max
h_rotor_iron = h_max
h_shaft_iron = h_max

rotor_inner  = occ.Circle(origin,r=r1).Face()
rotor_outer  = occ.Circle(origin,r=r2).Face()
sliding_inner  = occ.Circle(origin,r=r4).Face()
sliding_outer  = occ.Circle(origin,r=r6).Face()
stator_inner = occ.Circle(origin,r=r7).Face()
stator_outer = occ.Circle(origin,r=r8).Face()

rotor_inner.edges[0].name = "rotor_inner"
rotor_outer.edges[0].name = "rotor_outer"
stator_inner.edges[0].name = "stator_inner"
stator_outer.edges[0].name = "stator_outer"

rotor_iron = rotor_outer - rotor_inner

air_gap_stator = stator_inner - sliding_outer
air_gap = sliding_outer - sliding_inner
air_gap_rotor = sliding_inner - rotor_outer

stator_iron = stator_outer - stator_inner

for k in range(48):
    (stator_coil,stator_air) = drawStatorNut(k)

    stator_coil.faces.name = "coil" + str(k)
    stator_air.faces.name = "air"

    stator_iron -= stator_coil
    stator_iron -= stator_air

    domains.append(stator_coil)
    domains.append(stator_air)

for k in range(8):
    (magnet1,air_magnet1_1,air_magnet1_2) = drawMagnet1(k)
    (magnet2,air_magnet2_1,air_magnet2_2) = drawMagnet2(k)

    magnet1.faces.name = "magnet" + str(2*k+1)
    magnet1.faces.maxh = h_magnets
    magnet1.edges[0].name = "magnets_interface"
    magnet1.edges[1].name = "magnets_interface"
    magnet1.edges[2].name = "magnets_interface"
    magnet1.edges[3].name = "magnets_interface"

    magnet2.faces.name = "magnet" + str(2*k+2)
    magnet2.faces.maxh = h_magnets
    magnet2.edges[0].name = "magnets_interface"
    magnet2.edges[1].name = "magnets_interface"
    magnet2.edges[2].name = "magnets_interface"
    magnet2.edges[3].name = "magnets_interface"

    air_magnet1_1.faces.name = "rotor_air"
    air_magnet1_2.faces.name = "rotor_air"
    air_magnet2_1.faces.name = "rotor_air"
    air_magnet2_2.faces.name = "rotor_air"

    air_magnet1_1.faces.maxh = h_air_magnets
    air_magnet1_2.faces.maxh = h_air_magnets
    air_magnet2_1.faces.maxh = h_air_magnets
    air_magnet2_2.faces.maxh = h_air_magnets

    rotor_iron -= magnet1
    rotor_iron -= air_magnet1_1
    rotor_iron -= air_magnet1_2
    rotor_iron -= magnet2
    rotor_iron -= air_magnet2_1
    rotor_iron -= air_magnet2_2

    domains.append(magnet1)
    domains.append(magnet2)
    domains.append(air_magnet1_1)
    domains.append(air_magnet1_2)
    domains.append(air_magnet2_1)
    domains.append(air_magnet2_2)

stator_iron.faces.name = "stator_iron"
stator_iron.faces.maxh = h_stator_iron

air_gap_stator.faces.name = "air_gap_stator"
air_gap_stator.faces.maxh = h_air_gap

air_gap.faces.name = "air_gap"
air_gap.faces.maxh = h_air_gap

air_gap_rotor.faces.name = "air_gap_rotor"
air_gap_rotor.faces.maxh = 1.08*h_air_gap

rotor_iron.faces.name = "rotor_iron"
rotor_iron.faces.maxh = h_rotor_iron

shaft_iron = rotor_inner
shaft_iron.faces.name = "shaft_iron"
shaft_iron.faces.maxh = h_shaft_iron

domains.append(shaft_iron)
domains.append(rotor_iron)
domains.append(air_gap_stator)
domains.append(air_gap)
domains.append(air_gap_rotor)
domains.append(stator_iron)

geo = occ.Glue(domains)

if pizza_slice:
    pizza = occ.MoveTo(*origin).Line(1).Rotate(90).Line(1).Close().Face()

    geo = pizza*occ.Glue(domains)

    geo.edges[0].name = "left"
    geo.edges[4].name = "left"
    geo.edges[24].name = "left"
    geo.edges[28].name = "left"
    geo.edges[32].name = "left"
    geo.edges[96].name = "left"

    geo.edges[2].name = "right"
    geo.edges[6].name = "right"
    geo.edges[26].name = "right"
    geo.edges[30].name = "right"
    geo.edges[46].name = "right"
    geo.edges[98].name = "right"

    rot = occ.Rotation(occ.Axis((0,0,0), occ.Z), 45)
    geo.edges[2].Identify(geo.edges[0], "per", 0, rot)
    geo.edges[6].Identify(geo.edges[4], "per", 0, rot)
    geo.edges[26].Identify(geo.edges[24], "per", 0, rot)
    geo.edges[30].Identify(geo.edges[28], "per", 0, rot)
    geo.edges[46].Identify(geo.edges[32], "per", 0, rot)
    geo.edges[98].Identify(geo.edges[96], "dsa", 0, rot)

geoOCC = occ.OCCGeometry(geo, dim=2)
geoOCCmesh = geoOCC.GenerateMesh()

import ngsolve as ng
ngsolvemesh = ng.Mesh(geoOCCmesh)
ngsolvemesh.Refine()
# ngsolvemesh.Refine()
# geoOCCmesh.SecondOrder()

from netgen.webgui import Draw as DrawGeo
DrawGeo(geo)

BaseWebGuiScene

#ngsolvemesh.Refine()

mesh = ngsolvemesh
plist = []
for pair in mesh.ngmesh.GetIdentifications():
    plist += list(mesh.vertices[pair[0]-1].point) + [0]
    plist += list(mesh.vertices[pair[1]-1].point) + [0]

from ngsolve.webgui import Draw as DrawMesh
DrawMesh(mesh, objects=[{"type" : "lines", "position" : plist, "name": "identification", "color": "purple"}])

BaseWebGuiScene

#####################################################################################################################
Mperp_mag1 = np.array([-0.507223091788922, 0.861814791678634])
Mperp_mag2 = np.array([-0.250741225095427, 0.968054150364350])
Mperp_mag3 = (-1)*np.array([-0.968055971101187, 0.250734195544481])
Mperp_mag4 = (-1)*np.array([-0.861818474866413, 0.507216833690415])
Mperp_mag5 = np.array([-0.861814791678634, -0.507223091788922])
Mperp_mag6 = np.array([-0.968054150364350, -0.250741225095427])
Mperp_mag7 = (-1)*np.array([-0.250734195544481, -0.968055971101187])
Mperp_mag8 = (-1)*np.array([-0.507216833690415, -0.861818474866413])
Mperp_mag9 = np.array([0.507223091788922, -0.861814791678634])
Mperp_mag10 = np.array([0.250741225095427, -0.968054150364350])
Mperp_mag11 = (-1)*np.array([0.968055971101187, -0.250734195544481])
Mperp_mag12 = (-1)*np.array([0.861818474866413, -0.507216833690415])
Mperp_mag13 = np.array([0.861814791678634, 0.507223091788922])
Mperp_mag14 = np.array([0.968054150364350, 0.250741225095427])
Mperp_mag15 = (-1)*np.array([0.250734195544481, 0.968055971101187])
Mperp_mag16 = (-1)*np.array([0.507216833690415, 0.861818474866413])

Mperp_mag = np.c_[Mperp_mag1,Mperp_mag2, Mperp_mag3, Mperp_mag4, Mperp_mag5, Mperp_mag6, Mperp_mag7, Mperp_mag8,
                  Mperp_mag9,Mperp_mag10,Mperp_mag11,Mperp_mag12,Mperp_mag13,Mperp_mag14,Mperp_mag15,Mperp_mag16]

m = np.c_[Mperp_mag[1,:],-Mperp_mag[0,:]].T
nu0 = 10**7/(4*np.pi)
m = m*nu0*1.158095238095238
#####################################################################################################################




#####################################################################################################################
offset = 0
polepairs  = 4
gamma_correction_model = -30.0
gamma = 40.0
gamma_correction_timestep = -1
phi0 = (gamma + gamma_correction_model + gamma_correction_timestep * polepairs) * np.pi/180.0



def f48(s):
    return (s-1)%48

area_coils_UPlus  = np.r_[f48(offset+1) , f48(offset+2)]
area_coils_VMinus = np.r_[f48(offset+3) , f48(offset+4)]
area_coils_WPlus  = np.r_[f48(offset+5) , f48(offset+6)]
area_coils_UMinus = np.r_[f48(offset+7) , f48(offset+8)]
area_coils_VPlus  = np.r_[f48(offset+9) , f48(offset+10)]
area_coils_WMinus = np.r_[f48(offset+11), f48(offset+12)]

for k in range(1,4):
    area_coils_UPlus = np.r_[area_coils_UPlus, f48(k*12+offset+1)]
    area_coils_UPlus = np.r_[area_coils_UPlus, f48(k*12+offset+2)]

    area_coils_VMinus = np.r_[area_coils_VMinus, f48(k*12+offset+3)]
    area_coils_VMinus = np.r_[area_coils_VMinus, f48(k*12+offset+4)]

    area_coils_WPlus = np.r_[area_coils_WPlus, f48(k*12+offset+5)]
    area_coils_WPlus = np.r_[area_coils_WPlus, f48(k*12+offset+6)]

    area_coils_UMinus = np.r_[area_coils_UMinus, f48(k*12+offset+7)]
    area_coils_UMinus = np.r_[area_coils_UMinus, f48(k*12+offset+8)]

    area_coils_VPlus = np.r_[area_coils_VPlus, f48(k*12+offset+9)]
    area_coils_VPlus = np.r_[area_coils_VPlus, f48(k*12+offset+10)]

    area_coils_WMinus = np.r_[area_coils_WMinus, f48(k*12+offset+11)]
    area_coils_WMinus = np.r_[area_coils_WMinus, f48(k*12+offset+12)]


I0peak = 1555.63491861 ### *1.5
phase_shift_I1 = 0.0
phase_shift_I2 = 2/3*np.pi#4/3*pi
phase_shift_I3 = 4/3*np.pi#2/3*pi

I1c = I0peak * np.sin(phi0 + phase_shift_I1)
I2c = (-1)* I0peak * np.sin(phi0 + phase_shift_I2)
I3c = I0peak * np.sin(phi0 + phase_shift_I3)

areaOfOneCoil = 0.00018053718538758062

UPlus  =  I1c* 2.75 / areaOfOneCoil
VMinus =  I2c* 2.75 / areaOfOneCoil

WPlus  =  I3c* 2.75 / areaOfOneCoil
UMinus = -I1c* 2.75 / areaOfOneCoil

VPlus  = -I2c* 2.75 / areaOfOneCoil
WMinus = -I3c* 2.75 / areaOfOneCoil

j3 = np.zeros(48)
j3[area_coils_UPlus]  = UPlus
j3[area_coils_VMinus] = VMinus
j3[area_coils_WPlus]  = WPlus
j3[area_coils_UMinus] = UMinus
j3[area_coils_VPlus]  = VPlus
j3[area_coils_WMinus] = WMinus
#####################################################################################################################

import sys
sys.path.insert(0,'../../../') # adds parent directory

import pde
import plotly.io as pio
pio.renderers.default = "notebook"


MESH = pde.mesh.netgen(geoOCCmesh)

ind_air_all = np.flatnonzero(np.core.defchararray.find(MESH.regions_2d,'air')!=-1)
ind_stator_rotor = np.flatnonzero(np.core.defchararray.find(MESH.regions_2d,'iron')!=-1)
ind_magnet = np.flatnonzero(np.core.defchararray.find(MESH.regions_2d,'magnet')!=-1)
ind_coil = np.flatnonzero(np.core.defchararray.find(MESH.regions_2d,'coil')!=-1)
ind_shaft = np.flatnonzero(np.core.defchararray.find(MESH.regions_2d,'shaft')!=-1)

trig_air_all = np.where(np.isin(MESH.t[:,-1],ind_air_all))
trig_stator_rotor = np.where(np.isin(MESH.t[:,-1],ind_stator_rotor))
trig_magnet = np.where(np.isin(MESH.t[:,-1],ind_magnet))
trig_coil = np.where(np.isin(MESH.t[:,-1],ind_coil))
trig_shaft = np.where(np.isin(MESH.t[:,-1],ind_shaft))

vek = np.zeros(MESH.nt)
vek[trig_air_all] = 1
vek[trig_magnet] = 2
vek[trig_coil] = 3
vek[trig_stator_rotor] = 4
vek[trig_shaft] = 3.6

fig = MESH.pdemesh()
fig = MESH.pdesurf(vek)
fig.show()

def makeIdentifications(MESH):

    a = np.array(MESH.geoOCCmesh.GetIdentifications())

    femsol = np.zeros(MESH.np)
    femsol[a[:,0]-1] = +1
    femsol[a[:,1]-1] = -1

    fig = MESH.pdesurf(femsol)
    fig.show()

    c0 = np.zeros(a.shape[0])
    c1 = np.zeros(a.shape[0])

    for i in range(a.shape[0]):
        point0 = MESH.p[a[i,0]-1]
        point1 = MESH.p[a[i,1]-1]

        c0[i] = point0[0]**2+point0[1]**2
        c1[i] = point1[0]**2+point1[1]**2

    ind0 = np.argsort(c0)

    aa = np.c_[a[ind0[:-1],0]-1,
               a[ind0[1: ],0]-1]

    edges0 = np.c_[a[ind0[:-1],0]-1,
                   a[ind0[1: ],0]-1]
    edges1 = np.c_[a[ind0[:-1],1]-1,
                   a[ind0[1: ],1]-1]

    edges0 = np.sort(edges0)
    edges1 = np.sort(edges1)

    edgecoord0 = np.zeros(edges0.shape[0],dtype=int)
    edgecoord1 = np.zeros(edges1.shape[0],dtype=int)

    for i in range(edges0.shape[0]):
        edgecoord0[i] = np.where(np.all(MESH.EdgesToVertices[:,:2]==edges0[i,:],axis=1))[0][0]
        edgecoord1[i] = np.where(np.all(MESH.EdgesToVertices[:,:2]==edges1[i,:],axis=1))[0][0]

    identification = np.c_[np.r_[a[ind0,0]-1,MESH.np + edgecoord0],
                           np.r_[a[ind0,1]-1,MESH.np + edgecoord1]]
    ident_points = np.c_[a[ind0,0]-1,
                         a[ind0,1]-1]
    ident_edges = np.c_[edgecoord0,
                        edgecoord1]
    return ident_points, ident_edges

ident_points,ident_edges = makeIdentifications(MESH)
print(ident_points,ident_edges)

[[    1     1]
 [ 3360  8310]
 [   79    77]
 [12245  8309]
 [    2     0]
 [ 3363  3365]
 [  383    91]
 [12273  8324]
 [  382    90]
 [ 3781  8834]
 [  381    89]
 [ 3779  3376]
 [  380    88]
 [12272  3375]
 [  379    87]
 [ 3776  3849]
 [  378    86]
 [ 8771  3373]
 [  377    85]
 [ 8770  3372]
 [  376    84]
 [12279  8833]
 [  375    83]
 [12270  8320]
 [  374    82]
 [12269  8317]
 [  373    81]
 [12278  8316]
 [  372    80]
 [12251  8315]
 [    4     3]
 [10024 10023]
 [   22    21]
 [10472 10470]
 [   24    23]
 [10630 13078]
 [   38    25]
 [11261 11372]
 [ 1450  1434]
 [11371  6961]
 [ 1449  1433]
 [ 6880 11363]
 [ 1448  1432]
 [ 6879 11405]
 [ 1447  1431]
 [11591 11362]
 [ 1446  1430]
 [ 7335 11361]
 [ 1445  1429]
 [13175 13174]
 [ 1444  1428]
 [ 6878  6867]
 [ 1443  1427]
 [ 6877  6866]
 [ 1442  1426]
 [11667  6865]
 [   76    75]] [[    5     6]
 [  584   574]
 [  587   573]
 [   11     3]
 [    9     2]
 [ 2921   642]
 [ 2925   644]
 [ 2920   639]
 [ 2915   640]
 [ 2912   634]
 [ 2911   631]
 [ 2908   628]
 [ 2910   627]
 [ 2906   623]
 [ 2902   624]
 [ 2898   619]
 [ 2900   618]
 [ 2895   614]
 [ 2894   613]
 [ 2890   609]
 [ 2891   612]
 [ 2887   608]
 [ 2886   607]
 [ 2882   601]
 [ 2881   600]
 [ 2876   596]
 [ 2877   595]
 [ 2872   591]
 [ 2871   590]
 [   21    14]
 [   20    16]
 [  154   149]
 [  156   150]
 [  164   159]
 [  165   161]
 [  270   171]
 [  272   170]
 [11326 11251]
 [11327 11249]
 [11323 11244]
 [11319 11246]
 [11316 11241]
 [11315 11242]
 [11310 11237]
 [11314 11236]
 [11308 11231]
 [11307 11230]
 [11302 11226]
 [11304 11227]
 [11299 11222]
 [11297 11219]
 [11292 11214]
 [11291 11213]
 [11286 11208]
 [11289 11207]
 [  569   564]]

# import matplotlib.pyplot as plt
# plt.plot(np.sort(c0))
# MESH.pdemesh(info=1)

#####################################################################################################################
#import scipy.io
#scipy.io.savemat('motor.mat', {"m" : m,
#                               "j3" : j3,
#                               "geoOCC" : geoOCC,
#                 do_compression=True)
#####################################################################################################################


if pizza_slice:
    #####################################################################################################################
    np.savez_compressed('motor_pizza.npz', m=m, j3=j3, geoOCC = geoOCC,
                        geoOCCmesh = geoOCCmesh, ident_points = ident_points, ident_edges = ident_edges)
    #####################################################################################################################
else:
    #####################################################################################################################
    np.savez_compressed('motor.npz', m=m, j3=j3, geoOCC = geoOCC,
                        geoOCCmesh = geoOCCmesh)
    #####################################################################################################################

ident_points

array([[    1,     1],
       [ 3360,  8310],
       [   79,    77],
       [12245,  8309],
       [    2,     0],
       [ 3363,  3365],
       [  383,    91],
       [12273,  8324],
       [  382,    90],
       [ 3781,  8834],
       [  381,    89],
       [ 3779,  3376],
       [  380,    88],
       [12272,  3375],
       [  379,    87],
       [ 3776,  3849],
       [  378,    86],
       [ 8771,  3373],
       [  377,    85],
       [ 8770,  3372],
       [  376,    84],
       [12279,  8833],
       [  375,    83],
       [12270,  8320],
       [  374,    82],
       [12269,  8317],
       [  373,    81],
       [12278,  8316],
       [  372,    80],
       [12251,  8315],
       [    4,     3],
       [10024, 10023],
       [   22,    21],
       [10472, 10470],
       [   24,    23],
       [10630, 13078],
       [   38,    25],
       [11261, 11372],
       [ 1450,  1434],
       [11371,  6961],
       [ 1449,  1433],
       [ 6880, 11363],
       [ 1448,  1432],
       [ 6879, 11405],
       [ 1447,  1431],
       [11591, 11362],
       [ 1446,  1430],
       [ 7335, 11361],
       [ 1445,  1429],
       [13175, 13174],
       [ 1444,  1428],
       [ 6878,  6867],
       [ 1443,  1427],
       [ 6877,  6866],
       [ 1442,  1426],
       [11667,  6865],
       [   76,    75]])

geoOCCmesh.Points()

<netgen.libngpy._meshing.Array_class netgen::MeshPoint_class netgen::PointIndex at 0x183dabeeaf0>